function epochNew=bciNoiseCancellation(dat,reference)
% epochsNew=noiseCancellation(data,reference,outFNname) MEG-drift removal for BTI-systems using reference channels 
% INPUT:
% dat:     Array with data to be correct
%           format: channel,time
% refernce: Array containing the reference channel signals
%           format: channel,time
%           Note that not all reference channels must be passed
%           A subset may even produce better results (e.g. we use the
%           reference channels 1-9 for our 248 channel system.
%           Note that time and trial inidces must correspond with "data"

% Written and developed by Markus Heiny and Jochem Rieger, Dept. of Neurology II, Otto-von-Guericke University Magdeburg
% For questions please contact: jochem.rieger@med.ovgu or markus@r-heiny.de
% Code based on Robinson, S.E. (1989) Noise Cancellation For Biomagnetic Measurements. In: Advances in Biomagnetism, Plenum Press, New York 1989 
% 20111115 CR: code adapted to work with single epochs only


%M=number of signalchannels
%N=number of reference channels
%L=length of the timevector

% W=zeros(size(reference,1),size(data,1));%                                                     NxM
s0=dat-repmat(mean(dat,2),1,size(dat,2));     %mean set to zero                                      MxL
r0=reference-repmat(mean(reference,2),1,size(reference,2));     %mean set to zero                                      NxL

rcov=r0*r0';                            %covariance of the reference channels                    NxN
rq=diag(rcov);                        %variance of the reference channels                    N

rs=r0*s0' ;                             %covariance between reference- and  signalchannels     NxM

 W=rs./repmat(rq,1,size(rs,2));          %weightmatrix                                          NxM
% W=rs./(rq*ones(1,size(rs,2))); %alternative calculation

 U=rcov./repmat(rq',size(r0,1),1);   %orthogonalizationsmatrix                              NxN
% U=rcov./(ones(size(r0,1),1)*rq'); %alternative calculation

epochNew=dat-W'/U*r0;

